﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.NonlinearEquationSolvers
{
    /// <summary>
    /// This class provides the finding of roots of a function by using the Brent-Dekker method.
    /// </summary>
    [Serializable]
    public class BrentMethodRootFinder : AbstractDerivativeNoNeedRootFinder
    {
        /// <summary>
        /// Initializes a new instance of the <see cref="BrentMethodRootFinder"/> class.
        /// </summary>
        /// <param name="polynomial">The polynomial for finding the roots.</param>
        public BrentMethodRootFinder(Polynomial polynomial)
            : base(polynomial)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="BrentMethodRootFinder"/> class.
        /// </summary>
        /// <param name="polynomial">The polynomial for finding the roots.</param>
        public BrentMethodRootFinder(SimplePolynomial polynomial)
            : base(polynomial)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="BrentMethodRootFinder"/> class.
        /// </summary>
        /// <param name="function">The function for finding the roots.</param>
        public BrentMethodRootFinder(IRealFunction function)
            : base(function)
        {
        }

        /// <summary>
        /// Finds one root of the polynomial by using Brents method. If the function
        /// has not real roots, which are between the specified values a and b, the result of the
        /// method is wrong. So the result of the method has to checked by the property PrecisionError.
        /// </summary>
        /// <param name="a">The x value for the first starting point. It has to be a value on the left
        /// side of the root to find.</param>
        /// <param name="b">The x value for the second starting point. It has to be a value on the right
        /// side of the root to find.</param>
        /// <returns>
        /// A root of the polynomial, which lies between the specified values a and b.
        /// </returns>
        public double FindRoots(double a, double b)
        {
            return this.FindRoots(a, b, 1e-15, 1000);
        }

        /// <summary>
        /// Finds one root of the polynomial by using Brents method. If the function
        /// has not real roots, which are between the specified values a and b, the result of the
        /// method is wrong. So the result of the method has to checked by the property PrecisionError.
        /// </summary>
        /// <param name="a">The x value for the first starting point. It has to be a value on the left
        /// side of the root to find.</param>
        /// <param name="b">The x value for the second starting point. It has to be a value on the right
        /// side of the root to find.</param>
        /// <param name="iterations">The maximum of iterations to do, for finding the root.</param>
        /// <returns>
        /// A root of the polynomial, which lies between the specified values a and b.
        /// </returns>
        public double FindRoots(double a, double b, int iterations)
        {
            return this.FindRoots(a, b, 1e-15, iterations);
        }

        /// <summary>
        /// Finds one root of the polynomial by using Brents method. If the function 
        /// has not real roots, which are between the specified values a and b, the result of the 
        /// method is wrong. So the result of the method has to checked by the property PrecisionError.
        /// </summary>
        /// <param name="a">The x value for the first starting point. It has to be a value on the left 
        /// side of the root to find.</param>
        /// <param name="b">The x value for the second starting point. It has to be a value on the right 
        /// side of the root to find.</param>
        /// <param name="precision">The precision of the result.</param>
        /// <param name="iterations">The maximum of iterations to do, for finding the root.</param>
        /// <returns>A root of the polynomial, which lies between the specified values a and b.</returns>
        public double FindRoots(double a, double b, double precision, int iterations)
        {
            int i = 0;
            double fa = this.Function.SolveAt(a);
            double fb = this.Function.SolveAt(b);
            double c = a;
            double fc = fa;
            double d = b - a;
            double e = d;

            if (fa * fb > 0)
            {
                throw new ArgumentException("The values at f(a) and f(b) have to have different signes.");
            }

            while (i < iterations)
            {
                if (fb * fc > 0)
                {
                    c = a;
                    fc = fa;
                    d = b - a;
                    e = d;
                }

                if (Math.Abs(fc) < Math.Abs(fb))
                {
                    a = b;
                    b = c;
                    c = a;
                    fa = fb;
                    fb = fc;
                    fc = fa;
                }

                double m = (c - b) / 2;
                double toleranz = 2 * precision * Math.Abs(b);

                if ((Math.Abs(m) > toleranz) && (Math.Abs(fb) > 0))
                {
                    if ((Math.Abs(e) < toleranz) || (Math.Abs(fa) <= Math.Abs(fb)))
                    {
                        d = m;
                        e = m;
                    }
                    else
                    {
                        double q;
                        double p;
                        double s = fb / fa;

                        if (a == c)
                        {
                            p = 2 * m * s;
                            q = 1 - s;
                        }
                        else
                        {
                            double r = fb / fc;

                            q = fa / fc;
                            p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
                            q = (q - 1) * (r - 1) * (s - 1);
                        }

                        if (p > 0)
                        {
                            q = -q;
                        }
                        else
                        {
                            p = -p;
                        }

                        s = e;
                        e = d;

                        if ((2 * p < 3 * m * q - Math.Abs(toleranz * q)) && (p < Math.Abs(s * q / 2)))
                        {
                            d = p / q;
                        }
                        else
                        {
                            d = m;
                            e = m;
                        }
                    }

                    a = b;
                    fa = fb;

                    if (Math.Abs(d) > toleranz)
                    {
                        b += d;
                    }
                    else
                    {
                        if (m > 0)
                        {
                            b += toleranz;
                        }
                        else
                        {
                            b -= toleranz;
                        }
                    }
                }
                else
                {
                    break;
                }

                fb = this.Function.SolveAt(b);
                i++;
            }

            this.NeededIterations = i;
            this.PrecisionError = i >= iterations;
            this.RelativeError = Math.Abs((c - b) / 2);

            return b;
        }
    }
}